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Abstract 

We explore experimental procedure? for comparing the capabilities of complex discrete < vent 
service systems. Instead of measuring system capability by analyzing or simulating the system 
with a constant rate of arriving work system capability is measured as the maximum rate of 
work arrival for which the system has a steady state. Hence, we seek the arrival rate which 
causes the system to be at full capacity. This rate is arguably the best indication of the service 
system’s capability. We treat both work-conserving and non- work-conserving service systems, 

Using traditional aiid ispcCnih^cd mcdsuic? ufayM-ciu peifui malice. 

Introduction 

As industrial engineers, applied probabilists, sirnulationists, and systems analysts, we are often 
called upon to evaluate systems which service input traffic and produce finished products. These 
sytems are sometimes traditional queues or networks of queues, but are often systems with queue-like 
characteristics which cannot accurately be modeled as traditional queueing systems. In practice and 
in the literature, this evaluation is traditionally based on excercising a model- of the service system 
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by subjecting it to a stream of input traffic and estimating or calculating some expected system 
performance measure. 

We feel that this typical experimental design is lacking, and that the shortcomings stem from 
the arbitrary choice of the distribution of the input process. Especially problematic are cases where 
the service system being modeled does not currently exist, where worst-case behavior is sought, 
or where we wish to evaluate the system in situations which are not accessible for data collection. 
Practical service system analysis is interesting only when the service system is in an environment 
where the workload is high relative to the system’s capability to serve. In all that follows, we are 
interested in finding the intensity of the input process that taxes the service system to the extremes 
of its capabilities, and using this intensity as a measure to compare systems. 

1 The General Service System Model 

The service systems considered all have the following features: 

1. a centralized, controlable, nonlattice process which generates tasks at a rate A per unit time; 

2. tasks are admitted upon generation and processed by the system; 

3. a completed task is ejected from the system; 

4. the system has the capability to process as many as // tasks per unit time. 

We will call such systems Discrete Event Service Systems (DESSs), see figure 1. We will allow 
the system to create, combine, destroy, or absorb tasks - the DESS need not be work conserving. 
A DESS must, however, be mixed or open, as we are interested in how much work we can inject 
into the system without overburdening it. We may measure the performance of the system using 
traditional queueing measures such as the number of tasks resident in the system or some production 
cost, or we may opt to analyze some measures which are particular to the application at hand. In 
this work, we will approach the work-conserving and nonwork-conserving systems seperately. 
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Input process 
with intensity 

Mi) 



Output process 
witli maximum 
intensity p 



Figure 1: A simple DESS. 

1.1 Motivating Example 

Recently a mo<M was constructed of all of the single-channel radio communications in a Marine 
Expeditionary Brigade (MEB). A MEB consists of approximately ’20,000 marines who perform in 
three subelements, the Ground Combat Element (GCE), the Air Combat Element (ACE), and the 
Co in hat S^rvirp Support El^m^n* (CSSE). Th 4 ? single-channel radio network of the MEB may consist, 
of several thousand radios operat ing on several hundred radio nets. The goal of the study was to 
allocate newly purchased antijamming radios to the different units in the MEB. 

The radio traffic was modeled using the Marine Corps’ version of a structured traffic model. 
There are sixty classes of communication task packages, each representing a different mission that 
the MEB executes in battle. Each task package, called a Broad Operational Subtask (DOST), was 
composed of several tasks which were called Message Exchange Occurrances (MEOs). Each BOST 
contained between five (5) and sixteen (1G) MEOs, and each MEO has a set of other MEOs in the 
BOST which must be completed before it can be initiated. 

The MEOs are completed by traiismittting a message from the specified sender to each of the 
specified receivers. If this communication cannot be accomplished using the specified doctrinal radio 
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net, the sender will attempt to reach the receiver using other nets, relays through other radios, or 
using messengers. The user may also elect to delay the MEO some short time before attempting it 
again. The radios in the system 

• experience failures; 

• become jammed by scenario-specified jammers; 

• hold queues of MEOs with different priorities; 

• are capable of moving into and out of range of other radios; 

• can be used in voice or digital mode; 

• exit and join different radio nets as they move around or attempt to reroute MEOs. 

The antijamming radio fails less frequently and is harder to jam than the old radio. However, it is 
much more time consuming to enter a net with the new radio than with the old one. Finally, if an 
old radio tries to contact a new radio, the new radio must transit into a less capable mode in order 
to receive the MEO, then reenter his regular net of new radios. 

Our model features structured traffic being presented to a set of servicing capabilities which act 
semi-autonomously. The routing of the MEOs, sometimes using several transmissions to accomplish a 
single MEO, makes this model extremely difficult to analyze. To expedite the analysis and to acheive 
maximum flexibility and sponsor acceptance, an object-oriented simulation model was constructed 
to allow analysts to build MEB radio net -structures and test their capabilities against on another. 

This brings us to the search for the right rate of generstion of BOSTs in the system. This 
generation rate is clearly dependent on the pace and nature of the battle being experienced by the 
MEB. After some initial searching through volumes of data from the recent Desert Storm operation, 
we found that the Marine Corps didn’t, take time during their ground war to record the time and type 
of every DOST they executed! Experience with the model showed that the measured performance 
depended greatly on the pace of the presented traffic. When considering various C 3 I architectures, 
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we found the ranking of the radio allocations from best to worst changed as the BOST generation 
rate changed. The sponsor wanted to identify the best architecture for the most intense traffic. 

2 Work Conserving Systems 

We first study the behavior of systems where there is a one-to-one correspondence between the 
tasks we submit for processing and the finished tasks the service system produces. Work-conserving 
queueing models do not allow 

• tasks to expire while in service; 

• tasks to create other tasks while in service; 

• tasks to be split or combined; 

• tasks which never finish service. 

Work-conserving queueing system models are common in both the practice and literature of 
applied probability. In a typical experiment, we generate input to the system at a constant rate, 
monitor the performance of the system either at fixed intervals or upon departure from the system, 
and employ well-known methods of steady-state analysis to estimate the steady-state average of the 
performance measure. 

A maxim of the analysis of service systems is that the system will have stationary long-run 
behavior if and only if the number of arriving tasks are, on average, less than the number of tasks 
the system is capable of processing. If our overall system can work at a maximum of // tasks per 
unit time, we can input as many as // per unit time and the system will remain stationary.. If A is 
our arrival rate for the system, we wish to manipulate A to expose //. 

2.1 Generating Data 

There are two ways we can generate data from a work-conserving system which will reveal the 
maximum processing rate in the system. They are: 
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• input, tasks to the system at a rate known to be much higher than the system can handle; 

• fill the system, then input a new task every time that a task completes. 

Instead of choosing a very high input rate and dealing with the problems of exploding buffer contents 
and a nonrecurrent system, we will simply close off the system and recirculate the tasks which finish. 
This approach will also serve as a good introduction to the analysis of systems which do not conserve 
work . 

Thus, we examine a special kind of closed queueing network - one with a single loop-back which 
all tasks traverse. Let A (t) be the time-dependent rate of recirculation of tasks in the system. So 
long as the system contains enough tasks to keep it working at capacity, we have A(/) — ► // as 
/ — oo. Kelly [3] and Walrand [9] both show this for exponentially distributed service, and Disney 
and Kiessler [2] make the extension to Jackson networks. The result can be extended in the obvious 
way by treating Phase-type distributions for service times, to produce the result we seek (A (/) — ► fi 
as t oo) for generally distributed service. 

Example 

We will demonstrate this method on the Jackson network shown in figure 2. 

2.2 Detecting Transition to Steady State 

Let us simulate the completion of the first, N customers serviced by the closed system for M inde- 
pendent replications. Let T{ j be the time between recirculation during the i^ replication. Thus, 
T t j , i — 1, 2, . . . , M is a set of iid samples. Let Tj = Y2iLi Ti t j/M be the average recirculation time 
process. We seek the index A T * such that ETi tJ = ET = fi for all j > A r *. Hence, we are in the 
setting of a traditional initial transient detection problem. 

There exist many ways to tackle this problem, including 

• cross-replication confidence intervals, [10] 



tests for significant drift; 
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Figure 2: A Jackson Network designed to have a maximum service rate of 0.5. The numbers in 
parentheses are the number of servers at cacti station, and the routing probabilities are shown on 
the workstation connections. All servers have unit service time, and all bufTers are infinite. The 
dashed line shows the recirculation route added to force the system to serve at the maximum rate. 
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• standardized time series (STS), [6]. 

In our pxnerienrps wp Ii^vp found n«pfn1 pcn^rnljv in a slightly modified version we have 

developed, which we call ratio STS (RSTS). 



2.3 Ratio STS 



The method of standardized time series (STS) [6], produces confidence intervals from autocorrelated, 
stationary data. This method was used in [7] to detect the existence of initialization bias in 
simulation output, and was sharpened to produce optimal tests when the functional form of the 
initialization bias is known. 

Suppose that we have M independent samples of n points each, with Y ZJ being the j th point in 
the i th independent sample. Let 

Yij = r l j2 Yi * (1) 

* = i 

for i = 1,2,..., M , and with Y J>0 = 0 for each i. The time series Si(k), k = 1,2 , ; s constructed 

for each independent replication i as 



Si(k) 



Yi t n - Yi'k for 0 < k < n 
0 k = 0 , n . 



( 2 ) 



Let <t be the variance of V : . If S t (k) is divided by Vy/n/k and scale the index k so that the result 
resides in the unit interval [0, 1], the resulting time series 7i(/),0 < / < 1 is known to approximate a 
Brownian bridge as n — * oo. This is the fundamental result of [6], and the theoretical basis of this 
sequencial procedure. 

Schruben shows that scaling and summing T 2 (f), 

n 

= <J\fn T,(kn) (3) 

k = 1 

results in a normal random variable with variance given by 

a 2 7i(?i 2 — 1 ) 



VAR(Ai) = 



12 



(4) 
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Note that, except for a factor of <r 2 , VAR(A{) is independent of the data, it relies only on the 
parameters of the experiment. Hence, for any integer d < M , 

A 

A x 



\d 



B- 



1 1 \/VAR{A t ) 

12 



)‘ 



- 1Z Y ",! 2 

cr~n(n- — 1 ) ^ ' ' 



(5) 



( 6 ) 



Tiie original STS used to detect initial transients used as a test statistic for stationarity of the 
mean response. If we form a ratio of and \h-d> we can eliminate the need to estimate cr 2 , 
forming 






d.M —d 



i A; 



(M d) 1 Ylt = M-d^i 



(?) 



This test statistic, which we call the RSTS test statistic, is easy to use in all of the applications 
where STS is applied. In particular, if we are interested in determining the onset of steady state, we 
can form the backward-moving sequences j, j — n — l,n — 2, . . ., 1 for each replication 2 , where 
j is formed from the subsequence Y t k — j , j -h 1, . . . , ?t, the port, on of the i th replication between 



j and n. Thus, we form the sequence of F-statistics 



F d,M-d{j) 









1 J 



( 8 ) 



If we assume that the system is in steady state when each of the A tn are collected, then we can detect 
the transition of the system into st.eadv state bv looking at the first index N* where Fj A r *) 
exceeds the critical value for an F random variable with identical degrees of freedom. This method 
is demonstrated in the following example. 



2.3.1 Example 

Continuing with the work-conserving system example, suppose that we 

• start the system with 25 tasks enqued at workstation 1 at time 0.0; 

• simulate N — 500 customer recirculations; 



• replicate M = 20 times. 
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sample 

Figure 3: Trajectory of the mean process 7) and the confidence intervals. 

Figure 3 shows the trajectory of 1) and the associated confidence interval process for the first 
140 recirculations. Clearly, by sample N* = 80 we have passed the criteria for being in steady state 
according to Welsh’s cross-replication confidence interval method. Furthermore, we can see that any 
detectable slope in the mean process is neglegible. When tested for our 20 independent samples, the 
drift of tested to be insignificant (7/ 0 : no drift, has p-value « 0.4). 

The mean time between recirculations in 1 00, 101 , . . . , 140 was T — 0.56, (confidence interval 
(0.55330,0.56691)), clearly not as fast as the // = 0.50 which we know to be the system’s capacity. 
Performing unweighted RSTS in the first 140 samples showed no transition to steady state detectable 
- the procedure seemed to be accurate enough to discern that the transition had not yet occured, 
see figure 4. 
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Figure 4: RSTS performed on the first 140 sample recirculation times, using 12 numerator degrees 
of freedom and 8 denominator degrees of freedom Conclusion: no transition to steady state was 



evident.. 
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0 100 200 300 400 500 



sample 

Figure 5: RSTS performed on the first 500 sample recirculation times. 

When we extend the length of the runs we consider to the full 500 samples, we see that RSTS 
was able to indicate a strong transition to steady state around the N* — 330 sample, see figure 
5. Averaging the samples collected in 350,351, . . .,500, we observe an overall average of T — 0.501 
(confidence interval (0.49816,0.50389)). 

RSTS clearly dominated the other traditional initial transient methods. In the case of the 
recirculating jobs, we clearly have a very gradual descent to the steady-state average. The detection 
method is not important to our overall theme, though we must issue a general caution: The choice 
of N* should be made very conservatively. 
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3 Non- Work-Conserving Systems 

In this section, we consider the case where work is not conserved in the system, and where the 
measure of performance is very general. This case, which is much more interesting and applicable 
than the work-conserving case, allows us to treat cases where the service system may share unusual 
characteristics with the real system, and where the measure of performance of the system may be 
dictated by the study sponsor. 

3.1 Motivating Example, Revisited 

In the motivating example, BOSTs were input to the system. Depending on the DOST type, the 
BOST might 

• splinter into several communications tasks, which may splinter further at a later time; 

• require partial of full reassembly at different points: 

• expire after is has reached a certain age. 

Thus, the system clearly doesn’t conserve work. 

The sponsor was interested in specifying 

• a mix of different BOST types which the input was comprised of; 

• a time allotment for each type of BOST, and a time when the BOST expires and is removed 
from the system; 

• a one-time cost, by BOST type, assessed when the time allotment expires with the BOST still 
in process; 

All of these requirements w re imposed because of the need to assess the system’s ability to handle 
communications as diverse as artillery targetting and mission execution traffic, medical evacuation 
requests, situation reports, intelligence traffic, logistics and administrative communications, and 
communications allowing the radio nets to counter radio jamming. 
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For this application, we constructed a penalty process p(t ), which superimposed lateness and 
expiration penalties from all of the traffic in the system, and accumulated these penalties in [0,t]. 
We compared C 3 1 architectures based on the rate-of-climb of this penalty when traffic was inserted 
into the system. The sponsor really wanted to know the answer to the following question: ” Which 
architecture has the best peak-load performance?” Since the model was an idealization of the real 
system, and since wartime traffic load data is not available, we were faced with the problem of 
determining what traffic rate represented peak loading to the system. 

3.1.1 The Workload Ramp 

Given that we do not know the service rate of the system, we can try to produce an estimate 
by modulating the intensity of the system’s input and observing the effect this has on the system 
performance. One might attempt this by stepping through some reasonable intensity values, or 
by doing some sort of iterative search. After considering several alternatives, we decided that a 
nonhomogenius input process with gradually increasing intensity might be appropriate. This idea 
is similar to testing a stereo system for its ability to play loud music - we gradually turn up the 
volume, listening for the point where the music begins to become distorted. 

The mechanics of generating the ramping workload process and calculating the likelihood ratio 
for a generated sample path are now presented. The requirement is to generate a nonhomogeneous 
Poisson process with jumppoints . . , xn from an intensity function \{t) given by 

{ A (0) + rt t> 0 

( 9 ) 

0 t < 0 

where A(0) > 0 is the initial intensity and r is the rate of climb of the intensity function. In this 
presentation sign(r) is not specified, and is a point of interest in future research. Let N (t) be the 
number of tasks arriving during [0 ,<]. From the above equation, we have 



E[N(t)) 


(10) 


/ [A(0) + rs]ds 
Jo 


(11) 
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— -L r / 2 / 9 

• \ - / • / - 



( 12 ) 



Thus, if xj , 2 * 2 , • • • is generated via (9), then a(ij , £ 2 , • • • » ?n) is a Poisson |)rocess with rate 1, 
[1]. The scheme for generating our ramping workload process is given as 



ro = 0 . 0 , j = 1 

while NOT DONE 

generate U ~ U[0, l] 

Tj — Tj_i - ln(U) 
ij - a~'(r ; ) = 7 A l 0H y* (0>3 + g^I 
end while 



Algorithm 1: The generation of the ramping intensity workload process. 



Let G(x) — 1 — G(x) be the complement of the joint distribution of x \,X 2 ,. ■ , i\v- Then 



^x,|x,_i(0 — - 1 ^ ^ I x — 1 ] 

_ f -[a(x,^ 1 -f-t)-a(x2-l)] 

_ e -[A(0)(x l _ 1 +t) + r/2(f 1 _ 1 +t) 2 ]-[A(0)(f 1 _ 1 ) + r/2(f l _ 1 ) 2 l > 



(13) 

(14) 

(15) 



Thus, the conditional density function (/) is given by 

9i. |*._,(0 = [A(0) + r(/+f I _ 1 )]e-f A l°)‘ +r / 2 l ,J+ ^- , « (16) 

= [A(0) + r(i + i,_i)]e _[A(0)+r(t/2+i: '- ,)lt . (17) 

This last presentation of ^.Ix.-jCO highlights the nature of the density function, with leading con- 
stant given as A(0) -f ri^, the rate at the time of the generation epoch, and the exponent given as 
—t times A(0) -f r(ih_i -f X{)/ 2. 



When appropriately calibrated, the ramping workload process will drive the DESS into regimes in 
which it is underutilized, progressing to the point of total utilization, and then becoming saturated. 
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3.2 Detecting the Transition to Overloaded for a DESS 

Let our system performance measure for input intensity A and time t have expected value 

Our only assumptions about are that it is responsive to changes in A. it grows at a rate similar to 

a degree s polynomial when A < p, and it grows faster when A > p. 

Hence, by estimating the s + I st derivative of the performance measure (s + 1 because we would 
like to deal with a mean-zero sequence), we can produce a sequence with 

• constant mean when A < //, 

• some drift when A > //. 

Let t\02> • • • , Li be n evenly spaced points in t ime, where \(ti) is believed to be less than the service 
capacity // of the DESS, and A(/ n ) is believed to be much greater than //. Our method for estimating 
DESS capacity will 

1. replicate the system performance M times using the workload ramp as the input process, 

collecting data ), tj ) for the i th replication; 

2. form the s -f \ si derivative data ,- 5+ ^(A(/j), /;), j — s -f 1 , s + 2, . . . , 7i, / = 1 , 2, . . . , 72 using 
sequencial differences; 

3. perform a transition detection to determine the point N* where tyj 5 + 1 ^(A (tj) y tj) no longer 
have constant mean 0. 

We will propose and evaluate three methods for detecting the transition in the performance 
measure: 

• modified A'-charts; 

• RSTS; 

• adaptive regression splines. 

The application of each of these three methods will be made more difficult by a common feature of 
data we have collected - although the mean performance becomes constant after several differencing 
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WORKSTATION 


P[DESTRUCTION] 


P[CREATION] 


INSERTION STATION 


left 


0.2 


0.0 


- 


top 


0.1 


0.4 


2 


bottom 


0.2 


0.1 


4 


right 


- 


0.3 


2 



Table 1: Destruction and Creation of Tasks within the DESS. Tasks that are destroyed are not 
considered completed. 

operations, the variance still grows. Hence, our model of the data from the system when unsat urated 
will be Y iyj ,i = 1 , 2, . . . , A/, j — 1 , 2 ,..., 7?, which are mutually independent, and where = 0 is 
constant and <jy is unknown and assumed to v ary. 

Example, Continued 

We modified the Jackson network described in section 2 so that it was no longer work-conserving, 
and so that it exhibited properties similar to the communications system described above. Table 1 
shows what can happen to tasks in the system when they complete service at each of the nodes. 

For this system, we still have a maximum input rate of A = 2, as workstation 1 is not interfered 
with, and still has two servers serving at unit speed. The tasks for the system are of three classes 
All have identical processing speeds at the workstations, but each class pays a different price for 
waiting in each workstation buffer. Finally, each task can stay in the DESS cost-free for some period 
of time. After this delay, the cost per unit time is accumulated based on the buffer the task resides 
in. Each task becomes costless once it leaves the system. All of the data on task classes is in table 
2 . 

The system’s measure of performance is the cost accumulated from the beginning of the simu- 
lation. We subjected this system to a ramped workload process which started with insertion rate 
A(0) = 1.0 and climbing at a rate r = 0.00666. Thus, the system capacity is reached at time 
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CLASS 


FREE TIME 


COST(LEFT) 


COST(TOP) 


COST(BOTTO.M) 


COST(RIGIIT) 


1 


2.0 


4 


1 


0 


1 


2 


0.0 


1 


2 


3 


4 


3 


1.0 


3 


2 


1 


0 



Table 2: Free Time and Holding Costs for Task Classes. 



150 = A r *. We continued each experiment to time 300. 

As we see from figure 6, we have the expected properties of constant mean and growing variance 
for the second derivative of the accumulated costs when t < 150, and something else occurring when 
t > 150. 

3.2.1 Control Charting 

The most straightforward methodology comes from the field of quality assurance, [5], and involves 
the use of a sequencial hypothesis test. We wish to know the first time that we can conclude that 
it is no longer plausible that the response mean is EY — 0. Our growing varaince causes us to take 
one of two actions: 

• use the cross-replication sample standard deviation to form a confidence interval; 

• model the growth of the standard deviation and use the model when setting control limits. 

Using the second derivative data from our example system, we can see that the former method 
is inclusive because of constant false alarms (the lower control limit makes scleral visits above the 
x-axis during the trajectory), while using the modeled standard deviation creates an envelope which 
the mean stays inside even when we know the system is exhibiting drift, see figure 8. 

o n o n cmn r . t> i , • r r> i , * 

H.OXO iui i-Mjtectioii oi oaturatioii 

Let us return to the construction of the sequencial RSTS methods. In this case, we have two 



alterations to make: 
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Figure 6: Trajectories of the Second Derivative of the Accumulated Cost Samples. I he cent er line is 
the mean of M = 20 independent replications, while the surrounding lines are the upper and lower 



normal confidence intervals. 



sample 2nd derlv 
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Figure 7: The sample standard deviation of the responses and a regression model of the growing 



standard deviation. 
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Figure 8: Control limit detection methods. 
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1. t he process must be reversed to detect transitions out of steady state; 

2. the growing variability of the data violate the assumptions under which T,(/) on [0, 1] was 
shown to be a Brownian bridge - the underpinning of the method. 

Let the sequences A j , j — 1,2 , ...,j for each replication ?, where is formed from the sub- 
sequence Y l a- k — 1,2, ..., 7 . Thus, we are moving through the output sequence in the forward 
direction (opposite the usual STS for initialization bias). 

The growing variability of the output must be attacked directly. The model which fits the output 
with A < // is a mean-zero model with linearly increasing standard deviation. The expected number 
of input tasks, a(f), i s a l so quadratically increasing and is intimately linked to the growth of the 
preformance measure and its variability. If we collect samples of the cost function V(A,/) such 
that there are a constant expected number of input tasks per sample, we can avoid the problem of 
growing variability. Our data will still be formed by taking sequencial differences, but the intervals 
of sampling will contract. Figure 9 shows empirically that sampling this way produces data with 
constant standard deviation in our example. Proving that this technique works in all cases is 
impossible because of the breadth of our generality here. 

Let us establish the time interval sequence /i,0,...,/ n su ch that we expect to inject exactly 
C tasks into the system between and L + i,? = l,2,...,/i — 1. Hence, given we can compute 
U+i = t x -F At using 



a(t{ -f At) — «(/,) = C . 



(18) 



Using 12, we produce the equation 



A(0)A/ + r -{2t,At + A/ 2 ) - C = 0, 



(19) 



leading to 



-K + A(0)) + s/r^if + 2r/,-A(0) + A(0) 2 + 2rC 



(20) 



r 



If we sample the DESS cost function, starting with some initial point t\ and with some arbitrary 
value C, we produce a sample with constant standard deviation. Taking the sequencial differences 
as above, we produce a sequence which has the properties required to perform RSTS. 
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Figure 9: Sample standard deviations from data collected using the modified sampling points 



method. 
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Figure 10: Trajectory for RSTS for the modified sampling points data collection. Result: clear 
transition at time 123.0 

In our example system, we used the modified sampling points to produce the F-statistic trajectory 
shown in figure 10. The RSTS method gave a clear signal that the system lost steady-state at time 
123.0, where A (/,) = 1.82 and p — 0.910. We also tried RSTS with evenly spaced points, and acheived 
approximately the same result. From this experience, we conclude that straightforward RSTS is 
fairly robust with respect to fluctuations or growth in variability when the changes are coordinated 
as they are in this experiment. Analytical investigation has shown that these fluctuations do not 
cancel in the RSTS test statistic. 
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3.2.3 Adaptive Regression Splines 

Adaptive regression splines, especially multivariate regression splines, are the focus of intense basic 
research, see [8]. In our application, the regression spline required is especially simple, as the model 
is of a single independent variable, and we seek a single knot in the regression spline at the point 
where the zero-mean model departs from the data. Using the methods in Larson [4], we can derive 



the location of this single knot analytically. 

The model used is stated as 

[ fio + 0i (t — A 7 *) T t t t < A " 

Y(t)=l (21) 

y fio T 02 ( t — A r * ) T c t / > A * , 

where N m is still our point where the regression model changes. If we further prescribe that 

00 = 0 , ( 22 ) 

01 = 0, (23) 

we will be fitting the model with zero mean to the left of the single knot. Let 

SSEl = Y y2 (0; _ (24) 

t<A r * 

SSE r = Y 0’(0 - M (25) 

t>iV* 

SSE = SSE l + SSE R . (26) 



The method involves the optimal location of N* to minimize SSE. The procedure provided in [4] 
can be simplified for our application into a single-pass examination of the means of the data, but is 
valid only when a linear model is appropriate for the data to the right of N* . This method is valid 
when there is growing variability, as seen in our application. 

In our example, we fit the model in (21) with three parameters, then restricted it to the mean- 
zero model prescribed in (22-23). Figure 11 shows the two regression models. The three parameter 
model located a knot at time N* = 136.25, when the input rate is 1.91 and the traffic intensity is 
0.95. With* the one parameter restriction, the adaptive spline located the optimal knot at 128. 0. We 
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Figure 11: Adaptive regression spline wit h a single knot, three parameter and one parameter models. 



should note that the linear model does fit. fair I> well to the right of AT*, as one might expect of any 
queuing application with less capacity than required. 



4 Conclusion 



In this work, we have described the various problems which arise when we are attempting to de- 
termine the service capacity of a black box service system. We divided this investigation into two 
distinct parts, one dealing with simple queueing systems which conserve work. In this case, we 



REFERENCES 



‘27 



showed the advantages of closing the system so that output from the system was recirculated, as the 
time between recirculations converges to the system's service rate. We investigated ways to detect 
this convergence, and showed how difficult, this is in a simple Jackson network example. 

The second part of the exploration dealt with systems which 

• do not conserve work; 

• have their performance measured using a general holding cost mechanism or some other per- 
formance measure. 

In this case, we showed that if we modulated the input process using a ramped-intensity workload 
process, we could drive the system from underutilized to saturated. We explored three methods 
which unveil the point where this transition takes place, and demonstrated each on an example. 

The wider Significance of this work is the beginning of an exploration for empirical methods 
for determining the capacity of a service system. This exploration is done not by representing the 
system using a queuing model which we know how to analyze a prion , but by using a realistic model 
of the system and measuring its performance in terms the user has in mind. 
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